function dS = CRTBP_n_Thrust(t,S,u,T,c)
x = S(1);
y = S(2);
z = S(3);
vx = S(4);
vy = S(5);
% vz = S(6);
m = S(7);

r1 = ((x+u)^2+y^2+z^2)^(1/2);
r2 = ((x+u-1)^2+y^2+z^2)^(1/2);

n = length(S);
dS = zeros(n,1);

dS(1:3) = S(4:6);

dS(4) = x  -  (1-u)*(x+u)/r1^3  -  u*(x+u-1)/r2^3 + 2*vy;
dS(5) = y  -  (1-u)*y/r1^3 - u*y/r2^3 - 2*vx;
dS(6) = -(1-u)*z/r1^3-u*z/r2^3;
direc = cross([x,y,z],[0,0,-1])/norm(cross([x,y,z],[0,0,-1]));
dS(4:6) = dS(4:6) + T/m * direc';

dS(7) = -T/c;
end